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c/3 ' Abstract. We present a comprehensive study of the lattice restricted primitive model, 

i.e., a lattice gas consisting of an equal number of positively and negatively charged par- 
ticles interacting via on-site exclusion and a 1/r potential. On the cubic lattice, Monte 
Carlo simulations show a line of Ncel points separating a disordered, high-temperature 
n^ ' phase from a phase with global antiferromagnetic order. At low temperatures the 

rj . (high-density) ordered phase coexists with the (low-density) disordered phase. The 

Q ' Neel line meets the coexistence curve at a tricritical point, Tt ~ 0.14, pt ~ 0.4. A 

, ^, , simple mean-field analysis is in qualitative agreement with simulations. 

T-H ! 

^ : INTRODUCTION 

^ ' It has long been realized that the presence of an appreciable concentration of free 

Q . ions in a system of overall charge neutrality gives rise to a variety of thermodynamic 

G^ I features not found in un-ionized systems. For example, in the restricted primitive 

model (RPM), a system of charged hard-sphere anions and cations, all with equal 
c^ ■ charge magnitude |g| and diameter a, one has the famous limiting law established 

by the work of Debye and Hiickel [1] in 1923, 



TJ ■ Aex -p3 

C : . — ^ as p ^ 0, T fixed. (1) 

Here A"^^ is the Helmholtz free energy (excess to the ideal-gas condition), V 
'k>( I is the system volume, ks Boltzmann's constant, T the temperature, p the to- 

;h ' tal number density of anions and cations, and Tjj is the inverse Debye length, 

Vjja^ = Anpq^ /ksT. This is in marked contrast to an un-ionized fluid, for which 
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A'^^ /kBTV approaches a second virial-coefficient term proportional to p^ rather 
than p^/^. 

Resuhs pertaining to phase transitions and criticahty in the RPM did not emerge 
until much after the work of Debye and Hiickel. It was not until 1976 that a 
systematic study by Stell, Wu, and Larsen appeared with the conclusion that the 
RPM should exhibit a liquid-gas coexistence curve with a critical point [2]. They 
further concluded that the critical density is much lower than for a simple argon-like 
fluid model, but found that that the shape and location of the coexistence curve 
was very sensitive to the details of the approximations that must inevitably be 
used in such a theoretical study. Monte Carlo simulations have yielded results fully 
consistent with their conclusions, but only in the last few years have the studies 
of different groups converged toward common values for the critical density and 
temperature of the RPM [3-5]. 

As one of us has discussed in earlier studies [6,7], the RPM can be thought of as 
a spin system with a long-ranged antiferromagnetic interaction 

J{r) = -^, (2) 

re 

and one-dimensional spins Sj with |sj| = 1; 

(l)ijir) = -J{r)si-Sj . (3) 

(Since the "spins" are charges, we should perhaps refer to the interaction as anti- 
ferroelectric rather than antiferromagnetic.) The lattice-gas version of such a spin 
fluid is simply the spin-1 Ising model with |s| = 1, 0, or -1, with s = representing 
a vacant site. We shall refer to this model, which forms the subject of this report, 
as the lattice restricted primitive model (LRPM). 

For a lattice gas with a nearest-neighbor J(r) instead of a Coulombic J(r), the 
spin-1 Ising model is a Blume-Capel model [8,9], which becomes equivalent, at 
full occupancy, to a nearest-neighbor spin- ^/2 Ising model that is exactly-soluble 
on many two-dimensional lattices [10,11]. On the square lattice one finds a Neel 
point as one lowers the temperature in the absence of an external field. Moreover, 
although exact results are lacking in three dimensions, one expects for bipartite 
lattices such as the simple cubic, that as one lowers the density along the A-line 
of Neel points in the p-T plane, one encounters a tricritical point below which the 
lattice gas phase-separates into paramagnetic and antiferromagnetic phases. This 
immediately raises the question whether this remains true for a lattice gas with a 
Coulombic J(r). 

Several years ago, in order to better understand the LRPM, we initiated Monte 
Carlo simulations on a simple cubic lattice. Our results are fully consistent with 
the presence of a A-line of Neel points terminating in a tricritical point, below which 
(in T) two phases coexist. Preliminary results were reported in Ref. [12]; here we 
provide more extensive results as well as details of our simulation procedure. We 
have supplemented our simulations with a mean-field analysis that exploits the 



similarity between the LRPM and a nearest-neighbor lattice gas. H0ye and Stell 
[13] have also made a more general study of the spin-1 Ising model that included 
the solution of the mean-spherical approximation (MSA) for a range-parametrized 
J(r). As discussed by these authors, the MSA is not appropriate to study criticality 
and phase separation in the spin-1 model with Coulombic J(r) (although, as they 
show, it does yield the correct Debye-Hiickel limiting law). The H0ye-Stell study, 
however, goes beyond the MSA with a number of general observations that strongly 
support tricriticality for a Coulombic J(r). 

More recent studies made by Ciach and Stell [14] of ionic models that include 
the continuum-space and lattice RPM reveal that generically such models can be 
expected to manifest criticality, tricriticality, or both, depending upon the precise 
forn of their Hamiltonians. For example, the work suggests that some extended- 
core lattice models may well show both liquid-gas like criticality and, at higher 
densities, a A-line with associated tricriticality. For the RPM, the critical properties 
that follow from application of RG analysis are found to be in the Ising universality 
class. 



MODEL 

We consider a lattice gas of particles interacting via site exclusion (multiple 
occupancy forbidden) and a Coulomb interaction u{rij) = SiSj/vij, where Tj^ = 
|rj — Tjl is the distance separating the particles (located at lattice sites Fj and r^), 
and Sj = +1 or — 1 is the charge of particle i. Exactly half the A^ particles are 
positively charged; the remainder are negative. They are restricted to a simple 
cubic lattice oi V = L"^ sites, with periodic boundaries. We assume a lattice 
constant a of unity, and adopt units in which q'^/aks = 1,5' being the magnitude 
of the charge. In what follows we treat charge, length, energy and temperature as 
dimensionless. 

We note in passing that in lattice simulations of ionic systems an alternative 
definition of the potential is possible, i.e., u{rij) = SiSj^{rij), where $ is the 
Green's function for Poisson's equation on the lattice in question. While the simple 
1/r potential and the lattice Green's function differ somewhat at short distances, 
they have the same asymptotic (large r) behavior, and we expect (qualitatively) 
the same phase diagram in either case. 

In this study we restrict attention to the simple cubic lattice. This lattice, like 
the body-centered cubic, admits a decomposition into two sublattices (the sites on 
one sublattice having all nearest neighbors in the other sublattice), which obviously 
facilitates formation of an ordered state resembling an ionic crystal. Indeed, it was 
proven some time ago that the LRPM on the simple cubic lattice exhibits long- 
range order at sufficiently low temperatures and high fugacities [15]. The model 
presents, as we will show, a strong tendency to assume a NaCl-like ordered state at 
low temperatures. It would be interesting to study the model on the face-centered 
cubic lattice or another structure that frustrates antiferromagnetic ordering. 



MEAN-FIELD ANALYSIS 
The energy of the lattice restricted primitive model is 

u=\i:'^^ (4) 

where it is understood that the particles occupy distinct lattice sites. We consider 
the LRPM with independent variables T (temperature) and p = N/V. It is helpful 
to think of this system as a three-state (i.e., spin-1) antiferromagnetic Ising model 
with long-range interactions (J ~ I/''")- On a fully occupied lattice (p = 1) we 
characterize ordering by the sublattice charge disparity or staggered magnetization 
0. In the context of MFT, this allows us to replace Eq. (4) by the corresponding 
expression for a nearest-neighbor (NN) lattice gas, multiplied by a suitable factor. 
Consider first the disordered system. Since the neighbors of any given particle 
are equally likely to bear the same or opposite charges, the mean-field estimate for 
the energy is zero, just as for the NN system. Next, suppose there is sublattice 
ordering; let p^ be the fraction of sites in sublattice A occupied by positive particles, 
etc., and let 

(5) 

(6) 
so that 

^A _ ^A 

(7) 

and 

pt + P^ = pUp^^=p. (8) 

Using these, the mean-field estimate for the energy of a pair of nearest-neighbor 
sites is 

UNN = pIpI + ptp""^ - pip^ - p^pI 

= V0^. (9) 

Similarly, the MF estimate for the energy of a pair of second-neighbor sites (which 
must belong to the same sublattice) is +p^0^/\/2, and so on. Thus MFT gives the 
energy of the LRPM as 

^M.^^s;^. (10) 
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where S = +1(— 1) if sites i and j belong to the same (different) sublattices. The 
sum (including the factor of 1/2) is simply the electrostatic energy of an ionic 
crystal. For the simple cubic lattice, 

ii:H)! = -3K^, (u) 

2"^ Tii 6' ^ ^ 

as V ^ oo, where a = 1.7457 is the Madelung constant [16]. Thus the mean- 
field energy per site is Umf = — ttp^0^/2, while the corresponding expression for a 
system with nearest-neighbor interactions is — 3p^0^. In the context of MFT, then, 
we may replace the 1/r potential with a nearest-neighbor interaction. 

From the preceding analysis it appears that we may treat the LRPM, in MFT, as 
if it were a nearest-neighbor lattice gas, but with an energy (and temperature) scale 
that is smaller by a factor of a/6 ~ 0.2913. A more meaningful way of relating 
the temperature scales, however, is to compare the energy change attending an 
elementary excitation, in this case, the interchange of a nearest-neighbor pair in 
a perfectly ordered lattice. In the simple cubic lattice this raises the energy by 
^Unn = 20 (the nearest-neighbor interaction J = 1), since ten nearest-neighbor 
pairs have like charges after the exchange. For the ionic crystal (NaCl structure) the 
corresponding energy change is 4(q; — 1) ~ 2.9903, so that AUlrpm — 0.1495Af/Ar7v 
Since the mean-field critical temperature for the (fully occupied) nearest-neighbor 
lattice gas is 6 on the cubic lattice, the corresponding result for the LRPM is 
about 0.9. The best numerical estimate for the lattice gas is Tc = 4.5115; the 
corresponding LRPM value is 0.674. (Extrapolation of our simulation results to 
p = 1 yields T^ ^ 0.6.) 

We may now develop the MFT of the LRPM by studying the nearest-neighbor 
lattice gas, bearing in mind the difference in temperature scales explained above. 
To estimate the entropy as a function of p and 0, we note that the number of 
allowed configurations on a sublattice of V/2 sites is (V"/2)!/[A^+!A^_!A^„!] where iV+ 
(A^_) is the number of positive (negative) particles and N^ = [(y/2) — N^ — N_] 
the number of vacant sites on the sublattice. Using iV± = (1 ± (f))pV/4, etc., and 
Stirling's formula, we obtain the entropy per site. 



-plnp - (1-p) In(l-p) + pln2 - ^[(1 + 0) ln(l + 0) + (1 - 0) ln(l 
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Let f = u — Ts be the Helmholtz free energy per site. To investigate the possi- 
bility of a free energy minimum with nonzero sublattice ordering 0, we note that 
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A free energy minimum therefore implies that 
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lM.,„i±i. (14) 

T 1-0 ^ ' 

Expanding the r.h.s. as 20 + 20^/3 + - ■ ■, we see that solutions with nonzero exist 
for T < Tc = 6 and p > pc = T/Tc. For p ^ Pc 
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(15) 



Thus T = 6p is a hne of second-order phase transitions (a Neel or lambda line), 
separating the high-temperature phase {x = 0) from one with sublattice ordering. 
Next we investigate the stability of the uniform-density state. Defining fp = f / p 
as the free energy per particle, we may write the chemical potential as 

/i = /p + - , (16) 

P 

where the pressure p is given by 

p = p^m^\ =-Tln(l-p)-3pV. (17) 



(In evaluating the derivative we use Eq. (14) to eliminate {d(j)/dp)T-) The pressure 
in the disordered phase is that of a simple lattice gas; but for T < T^, {dp/dp)T 
suffers a discontinuity at pc- Using Eq. (15) we have that for p ^ Pc 



^C1 




T ^ /3p 



1-p Va 



9p --2 , (1^ 



so that 




T;p=pc 



nw:-2r (1^) 



Thus the inverse compressibility vanishes when pc = 1/3, which implies T = 2. This 
marks a tricritical point, at which the Neel line meets the coexistence curve. (For 
T < 2 the extension of the Neel line, T = 6p, is in fact a spinodal.) Recalling the 
temperature rescaling, our MET gives a tricitical temperature of 0.3 for the LRPM, 
or Tt = 0.225 if we use the best estimate {Tc = 4.5115) in place of the mean-field 
result {Tc = 6) for the simple lattice gas. Typical pressure-density curves are shown 
in Figure 1. 

We construct the coexistence curve by finding the densities pi and p2 for which 
p{pi,T) = p{p2,T) and p{pi,T) = p{p2,T). This is done numerically, through 
an iterative procedure. A first guess for p2 is the mean of p^ (the high-density 
spinodal, found by setting the r.h.s. of Eq.(18) to zero), and p^, the density at 
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FIGURE 1. Mean-field theory predictions for the pressure versus density at temperatures 0.449 
(upper curve), Tt = 0.299 (middle), and 0.224 (lower). 
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FIGURE 2. Mean-field prediction for the phase diagram of the LRPM. 



which the pressure equals p{pj), its value on the low-density spinodal. A first guess 
for pi is the density at which p{pi) = p{pm), given by the low-density isotherm, 
p = 1 — e~^/^. We then compare chemical potentials; if /x(p2) > /^(pi), our guess 
for p2 lies above the coexistence value, and vice-versa. By refining the estimate 
for p2, recalculating pi, and repeating the comparison, we rapidly converge on the 
coexisting densities. For temperatures T < 0.5 one may use = 1 — 2e^^'^^'^^'^ ~ 1 
to show that p2 ~ 1 — e~^^'^ and pi ~ 2e~^^'^^'^. The mean-field phase diagram for 
the LRPM (with temperature scaled such that Tt = 0.3) is shown in Figure 2. 

SIMULATION METHODS 

Electrostatic energy of a periodic array 

This subsection presents the multipole expansion as an efficient alternative to 
Ewald summation for evaluating the electrostatic energy of a periodic array of 
point charges. Our approach is closely related to the one employed by Ladd in 
simulations of ionic and dipolar fluids [17]. We require the electrostatic energy of 
an infinite, periodic array of cells containing point charges. Each cell harbors N/2 
positive charges and and equal number of negative charges; the position of particle 
i in cell j is 

Xij- = Xi + R^- , (20) 

where Xj is the position in the central cell, and 

Rj = a(n^x + Uj^y + n^z), (21) 

where n^ may assume any integer value. (The cell structure is that of a simple 
cubic lattice; the positions Xj may or may not be restricted to a lattice.) 

We write the electrostatic energy of the array as f/ = (1/2) J2iLi Qi^i, where $i 
is the potential at the position of charge i. When calculating <l>j we consider a cube 
of side a (oriented parallel to the periodic array), centered on charge i, and write 

$i = $,,0 + $^,a, (22) 

where <l>j_o is the contribution from the charges (other than i) within the cube, and 
$i,a is the potential at the center of the cube (i.e., the position of charge i) due 
to the infinite collection of periodic image cells. (These cells include the image of 
charge i, as they must, if they are to be neutral.) The lattice sum defining $j^a is 
conditionally convergent, and is often treated using the Ewald technique [18]. Here, 
however, we evaluate $j^a via a multipole expansion. The contribution to ^i^a due 
to cell j (j is a shorthand for the cell indices n^, Uy, and n^), is 

li2/ + i i?;+^ ^ ' 



Here Yi^ is the spherical harmonic, Rj = clJu"^ + n"^ + n^, cosOj = ariz/Rj, 
tan (f)j = Uy/rij:, and 

TV 

ql,m = J2^^'■'iYl*mi(^^^'Pi)^ (24) 

i=l 

is the (/, m) muhipole moment of (any) cell, with (rj, 6i, (pi) denoting the position of 
particle i relative to a convenient origin, for example the cell center. (The monopole 
moment goo vanishes by charge- neutrality.) 

Summing the contributions from all cells (except the central one) yields 

t^2/ + l^ rY' ^ ^ 

where the sum on j stands for a sum on Ux-, Uy, and n^ (not all zero). This sum over 
cells is independent of the charge distribution (which appears only in the qim), and 
so need only be evaluated once. Many terms vanish by symmetry. Thus all terms 
with / odd are zero since Yim is then an odd function of cos^, while the distribution 
of cells in the cubic lattice is symmetric about 6 = tt. Symmetry also implies that 
nonzero contributions have m = 4n {n integer), since for each cell with azimuthal 
angle 0, there are corresponding cells with angles + 7r/2, (p + n, and + 37r/2, and 
1 + e*™'^/^ + e*™'^ + e^*™'^/^ = 4 if m = An, and zero if m is an integer not divisible by 
4. (Angle is not defined at the poles, but Yim vanishes there anyway for m ^ 0.) 
Next we put the lattice sums that multiply the multipole moments into a conve- 
nient form. Using the relations 



Y, 



lm\ 



2l + l{l-mj: ^^,_^^ j^^ 



- Pr {cos 9)6'""", (26) 



\ An (/ + m)! 



(-) = (-i)™fq:^ Pn^)^ (27) 



and 



yi,_„(^,0) = (-i)™rz„(^,0), (2^ 

we can write Eq. (25) in the form. 



^ _ V L v- ^Kcos^,) , ^ {l-m)\ 

I even I j ^j m=4 \'' "•" '"')■ 






(29) 



where the sum on m is restricted to integer multiples of 4, with /' the largest such 
multiple < /, and "c.c." denotes the complex conjugate of the preceding term. The 
rescaled multipole moments are 

N 

ql,n. = T.9^riPricos9,)e-'^f^. (30) 

i=l 

Cubic lattice symmetry implies that the lattice sums are real, since the azimuthal 
angles associated with cells at a given R and 6 are either (0,7r/2, vr, 37r/2), or 
(7r/4, 37r/4, 57r/4, 77r/4), or take the form (±0, n/2 ± 0, tt ± 0, 3n/2 ± 0). Thus 
we may write 

/' 

$i,a = Y. Yl ^ImSlm, (31) 

I even m=4 

where /' is the largest multiple of 4 that is < /, and 

S»-E^^. (32) 



TDl+l 



and for m > 0, 



_ ^ (/ - m)! ^ Pl^{cos9j) cosm0j 
and the rescaled moments now take the form 

N 

li,ra = E qiriPricosOi) cosm0i. (34) 

If Sim is the sum for the integer lattice then Sim = Sim/a''^^- 

The expansion actually begins at / = 4. To see this, let {n^, Uy, n^j/j be the set 
of sites at distance R from the origin of the (integer) cubic lattice. Then 

S:^E^^-Ei-, S i(3cos^«-l). (35) 

j 3 R {nx,ny,nz}R 

The inner sum is 

_1_ Y: iH-nl-nl~nl) = 0, (36) 

{nx,ny,nz}ji 

since the cubic lattice is symmetric under permutations of x, y, and z. 

Applied to the NaCl lattice (positive and negative charges on the two sublattices), 
the expansion to order / = 10 yields an energy per particle of —aq'^/Aneoa with 



a = 1.747561. (The standard value is 1.747565 [16].) For this highly symmetric 
arrangement the / = 4 and 8 moments all vanish; the / = 10 contribution to the 
energy has a magnitude of about 0.03 that of the / = 6 term. The lattice sums Sim 
converge rather quickly; the following results were obtained by summing over all 
sites with R < 30: 

;S4o = 3.108219 ;S44 = 0.0370026 

Sqo = 0.5733293 ;S64 = -0.001592581 

^80 = 3.25929309 ^584 = 5.487025 x 10^^ ;S88 = 0.035665665 

Sio = 1.00922399 ^Sio 4 = -1.84839558 x 10"^ ^Sio 8 = -4.2786935 x 10"^ 

In the simulation routine we use these results to construct a lookup table giving 
the total contibution of a pair of charges at sites r and r' to the potential energy. 
The table entry includes the direct (minimum-image) contribution, the lattice-sum 
contribution obtained via multipole expansion, and a dipolar contribution to be 
explained in the following subsection. 



Dipole Correction Term 

In this very brief subsection we call attention to an absolutely crucial point in 
the simulation of systems (on or off-lattice) with slowly decaying forces. It is that 
the energy of the fully periodic array differs from the infinite lattice sum (taken, 
e.g., in spherical order, as above), by a term proportional to the square of the cell 
dipole moment: 

^=IT. ^^^^ = ^ E ^4<^.0 + <^^,a] + ^M' (37) 

where 

M = J2 Sir^ (38) 

i 

is the dipole moment of (any) cell. This result was derived and discussed in detail by 
de Leeuw, Perram, and Smith in 1980 [18]. It follows from a careful analysis of the 
lattice sum (summed in a particular order, using a convergence factor) applies to 
both the multipole respresentation of <I>i^a developed above, and to the more familiar 
Ewald expansion. We therefore add the term 27rSjSjr?/3L^ to the contribution of 
charges i and j to the potential energy. Some preliminary studies, which did not 
include this dipole term, yielded bizarre results, such as the tendency for a system 
with p ~ 1/4 to crystallize into a uniform-density body-centered cubic structure at 
low temperatures. 



Canonical Ensemble Simulations 

The majority of the data reported here were obtained in canonical ensemble sim- 
ulations following the usual Metropolis prescription, i.e., trial configurations that 
lower the energy are accepted with probability 1, while those attended by an energy 
increase Af/ > are accepted with probability e~^^l'^ . An initial configuration is 
generated by placing A^ particles onto the lattice at random, with double occupancy 
forbidden. This naturally yields a rather high-energy configuration that must be 
allowed to relax before the system attains equilibrium. We generally used a config- 
uration relaxed at some temperature T as the initial configuration for a study at 
a nearby temperature T' . The energy and order parameter (defined below) were 
monitored in order to check for relaxation. 

Trial configurations are generated by two kinds of "moves." The first is a single- 
particle move: a particle, chosen at random, is displaced by Ar = na;X + nj,y + n2Z, 
where n^^ Uy and n^ are independent, integer- valued random variables distributed 
uniformly on the set {— "m, ...,772}. (We generally used m = 2.) If the trial site is 
occupied, the move is rejected; otherwise the energy change is evaluated and the 
move accepted according to the Metropolis scheme. At low temperatures, each 
particle typically has one or more neighbors of the opposite charge; single-particle 
moves, which tend to disrupt these dipolar pairs, will have a low acceptance rate. 
We therefore alternated single-particle moves with pair moves. In this case we 
select a particle i at random, and check if one of the nearest-neighbor sites (also 
selected at random) harbors a particle. (Call it j; if the site is vacant the trial ends.) 
Given an occupied neighbor, particle i is displaced by Ar as above, and particle j 
is placed at a randomly chosen nearest-neighbor site of particle i. (Naturally, both 
trial sites must be vacant for the move to occur.) The move is accepted or rejected, 
as before, based on the Metropolis criterion. The results reported here represent 
averages over (typically) 2 to 10 of runs, each consisting of on the order of 2 x 10'^ 
lattice updates (i.e., moves per particle). A run of this sort requires 1-2 days of 
cpu time on a fast DEC alpha machine. 

In addition to the energy and the specific heat (obtained from the variance of 
the energy), we monitored the charge charge-charge correlation function, 

g,,{r)=AJ2{siO)si^)), (39) 

|x|=r 

and the density-density correlation function, 

^ppW = ^E(l^(o)lk(x)l). (40) 

|x|=r 

Here s(x) = 1 (-1) if site x is occupied by a positively (negatively) charged 
particle, and is zero otherwise. The normalization A is chosen such that Qpp -^ 1 
for r ^ 00. We also studied the sublattice order parameter, 0. Writing x = 
ix + jy + kz, we have 



= ^E((-ir''^"'M^'j>^))- (41) 

The order parameter is unity for a fully ordered (i.e., NaCl structure) lattice, and 
zero when the sublattices bear no net charge. 

In order to gauge the onset of phase separation, we divided the system into cells 
of 4 X 4 X 4 sites, and studied cell-occupancy histograms, P{p)- In an homogeneous 
system we expect P{p) to be unimodal, while a bimodal distribution indicates phase 
coexistence, and a broad or plateau-like distribution incipient phase separation. 



Grand Canonical Ensemble Simulations 

In an effort to obtain independent, and, one hopes, more reliable information on 
the coexistence curve, we performed simulations in the grand canonical ensemble. 
The procedure is complicated somewhat by the fact that we are obliged to maintain 
charge neutrality, so that a pair of particles (one positive, one negative) must be 
added or removed simultaneously. Let A = e^^^ be the pair fugacity, pp being the 
associated chemical potential. Then the grand partition function is 

5(/5, \V)= E A^/^ E ^-'""^'^ ^ (42) 

Areven {c}^ 

where the second sum is over the set of nonoverlapping configurations of N/2 pos- 
itive and N/2 negative particles on a lattice of V sites. 

Let C be a valid configuration for N particles, and C one for A^ -|- 2 particles, 
obtained by inserting particles (one positive, one negative) at some pair of vacant 
sites in C. (The insertion sites need not be nearest neighbors.) Let w{C^>C') be the 
transition rate, in a Monte Carlo simulation, for going from C to C (pair insertion), 
and w{C' ^>C) the rate for the reverse process (pair deletion). The detailed-balance 
condition, 

WJC^C) _ .-pyu{C')^U{C)\ /.ox 

w{C'^C) ' ^ ^ 

may be realized in various ways. In our simulations, pair addition and deletion steps 
are as follows. A fraction g/2 of the moves are addition attempts, and equal number 
are deletion attempts; the remainder are single-particle and pair displacements 
(no change in A^) as described in the preceding subsection. (In practice we used 
g = 0.9.) In an insertion attempt we choose two sites (anywhere in the system) at 
random. If either or both are occupied, the attempt fails; otherwise we tentatively 
place a positively charged particle at one site, and a negative particle at the other. 
Given the starting configuration C, the probability of realizing a particular trial 
configuration C is g/2V'^. The new configuration is accepted with probability 
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In a deletion attempt, we choose one of the positively charged particles at random, 
and similarly one of the negative particles, and tentatively delete them. Starting 
with configuration C (with A^ + 2 particles), the probability of realizing C as the 
trial configuration is 2g/{N + 2)^. In this case we accept the new configuration 
with probability 
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The transition rates w{C^C') = {g/2V^)Pins and w{C'^C) = [2g / {N + 2y]Pdei 
are readily seen to satisfy detailed balance. 

RESULTS 

In this section we summarize our simulation results for the phase diagram of the 
LRPM. We studied the energy, specific heat and order parameter as functions of 
temperature, for a series of 12 or so density values, ranging from about p = 0.02 to 
0.85. The data presented below were obtained in canonical ensemble simulations 
with system size L = 16 unless otherwise noted. 

Figure 3 shows the energy and Figure 4 the specific heat versus temperature 
for a number of different densities. All the specific heat curves show a peak at 
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FIGURE 3. Energy per particle versus temperature. Densities (left to right) p = 0.04, 0.08, 
0.16, 0.25, 0.35, 0.4, 0.5, 0.6, 0.75 and 0.85. 
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FIGURE 4. Specific heat versus temperature for densities as indicated. 

a temperature Ti that increases with density for small densities, and then seems 
to saturate at around T = 0.11 for p > 0.35. The height of this peak decreases 
rapidly with increasing density. For densities > 0.4, there is a second peak, at a 
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FIGURE 5. Loci of specific heat maxima (open and fiUed circles: L — 16; crosses: L = 20). 
The dashed hne denotes the position of the maximum of \d(f)/dT\. 
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FIGURE 6. Order parameter versus temperature for densities (left to right) 0.25, 0.35, 0.4, 0.5, 
0.6, 0.75, and 0.85. 



temperature T^ > Ti, which increases roughly hnearly with p. The locations of 
the specific heat peaks in the p-T plane are shown in Figure 5; these loci appear 
to delineate a Neel line and a coexistence curve. We examine their interpretation 
in what follows. 

It is tempting to associate the high-temperature specific heat peak with the onset 
of global sublattice ordering. Studies of the order parameter and of the charge- 
charge correlation function support this interpretation. In Figure 6 we plot the 
mean value of |0| versus temperature for various densities. In the thermodynamic 
limit, this quantity should be strictly zero in the disordered phase, and should in- 
crease continuously as T is reduced below T^. Here, of course, we must anticipate 
finite-size rounding of the critical singularity, if any exists. The data nevertheless 
give strong support for T/v marking an order- disorder transition. In fact the tem- 
perature at which \d(j)/dT\ takes its maximum value falls very close to that of the 
specific heat maximum, as shown in Figure 5. 

Further support is provided by the comparison between two system sizes (for the 
same density, p = 0.4) shown in Figure 7: the larger system presents a sharper 
variation in |0| in the vicinty of the apparent transition. The high-T specific heat 
peak, barely visible in the L = 16 data, is quite prominent for L = 20; both peaks 
shift to slightly higher temperatures for L = 20. All of this is consistent with the 
system exhibiting phase transitions at Ti and T/v- We shall therefore refer to the 
set of transitions commencing around p = 0.4 as the Neel line. 

Further insight into the phase transition at or near the high-T specific heat peak 
is afforded by the charge-charge correlation function Qqqij-). Figure 8 shows Qqqir) 
for three temperatures in the neighborhood of T/v, for p = 0.75. In all of the cases 



shown, Qqq exhlblts damped oscillations. At the highest temperature (T = 0.36) the 
oscillations evidently decay to zero; for T = 0.35 the decay appears to be slower, 
while for T = 0.34 they seem to persist. A plot of l^f^gl versus 1/r (Figure 9) yields 
a clearer idea of the large-r behavior; from this graph it is evident that Qqq decays 
to zero for T > 0.35 but not for T = 0.34. (The specific heat peaks at T = 0.325.) 
Thus the condition lifnr^oo\gqq{f^)\ > is a criterion for the ordered phase, and the 
limiting value of \gqq\ provides, in principle, an alternative definition of the order 
parameter. In practice, we are able to obtain with greater precision (smaller 
statistical uncertainty), and so retain it as a measure of global order. 

It is worth noting that in the neighborhood of the lambda line, the density- 
density correlation function Qpp shows little structure: it has a low-amplitude peak 
at unit distance and then decays rapidly to its asymptotic value (see Figure 10). 

We now turn to the line of low-temperature specific heat peaks. The mean-field 
theory results described above lead us to interpret this as a coexistence curve, so 
that Ti(p) marks a transition from a homogeneous system to coexisting high- and 
low-density phases, the former disordered, the latter with sublattice ordering. To 
better understand what is happening at Ti, we examine the occupancy histograms. 
Some typical results for p = 0.75 are shown in Figure 11. For T = 0.12 the distri- 
bution is unimodal, and peaked near p = 0.9. As the temperature is lowered, the 
distribution broadens and the peak moves to p = 1 and increases in amplitude, sig- 
naling the emergence of a high-density phase. For T = 0.11 a plateau has appeared 
for 0.2 < p < 0.7, while for T < 0.10 the distribution is bimodal. (Note also the 
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FIGURE 7. Upper graph: specific heat, p — 0.4. Filled circles: L — 16; open circles: L — 20. 
Lower graph: order parameter for p — 0.4; symbols as in upper graph. 
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FIGURE 8. Charge-charge radial distribution function, p = 0.75, for three temperatures as 
indicated. 




odd-even oscillations in P, especially prominent for T = 0.09. These presumably 
represent a tendency to local charge neutrality.) All through the transition region, 
the order parameter is very nearly unity; for this density, the transition near Ti 
is not associated with the onset of sublattice order. (Such order is already present 
in the uniform phase for Ti < T < T^.) 

Figure 12 shows a series of occupancy histograms for a much lower density, 
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FIGURE 9. The data of Figure 8 plotted versus 1/r. Lower set: T ^ 0.36; middle: T ^ 0.35; 
upper: T = 0.34. 




12 3 4 5 

r 

FIGURE 10. Density-density radial distribution function, p = 0.75, T = 0.34. 



p = 0.25. At the highest temperature shown, T = 0.12, the distribution peaks 
in the vicinity of p = 0.2 and decays rapidly for higher densities. For T = 0.115 the 
distribution is already considerably broader, and by T = 0.110 a plateau near p = 1 
has appeared. At T = 0.105 the distribution has become very broad, the main peak 
has shifted toward p = 0, and a secondary peak has appeared at p = 1, indicating 
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FIGURE 11. Cell-occupancy probabilities, p — 0.75. Upper graph: open circles: T ~ 0.12, 
filled circles T = 0.115; middle: open, T = 0.110, filled, T = 0.105; lower: open, T = 0.10, filled, 
T = 0.09. 
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FIGURE 12. Cell-occupancy probabilities, p — 0.25. Upper graph: open circles: T ~ 0.12, 
filled circles T = 0.115; middle: open, T = 0.110, filled, T = 0.105; lower graph: T = 0.10. 



separation into high- and low-density phases. By T = 0.10 the transformation is 
complete: the histogram peaks at p = and p = 1. 

It is interesting to observe the changes in the density-density correlation function 
hnn(r) = gpp{r) — 1, in the neighborhood of the transition. In Figure 13 we see that 



"ppy 



h decays rapidly for T = 0.12 (the correlation length ^p ^ 0.9), and is somewhat 
longer-ranged (^p ^ 2.2) for T = 0.115. For T = 0.11 the decay is much slower 
[^p ^6). At temperatures below 0.11, the system has phase-separated and we 
are no longer able to infer the large-r limiting value of gpp{r). The histogram and 
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FIGURE 13. Density-density correlation function, p ~ 0.25, T == 0.12 (lower set), 0.115 (mid- 
dle), and 0.11 (upper). 




FIGURE 14. A typical configuration, p ^ 0.12, T = 0.1. Filled and empty symbols distinguish 
the charges. 



correlation function data suggest that for p = 0.25, phase separation occurs at 
T = 0.105 — 0.110, in agreement with transition temperature of 0.108 indicated 
by the specific heat maximum. In contrast to what is seen at higher densities, for 
p = 0.25 the transition at Ti is accompanied by the onset of sublattice order: |0| is 
small for T > 0.12 and grows rapidly in the range T = 0.105 — 0.115. This is just 
what one would expect, since the transition in this case marks the appearance of a 
high-density phase, capable of supporting NaCl-like order. (Sublattice ordering is 
associated with the transition at Ti for p < 0.25.) 

Figure 14 is a snapshot of a typical configuration of the phase-separated system, 
in which we see a dense, crystal-like agglomerate coexisting with a dilute vapor, 
comprised largely of dipolar associations. 

Using the specific heat data, together with results for the order parameter and 
the charge-charge correlation function, we have derived a series of estimates for 
points along the lambda line, Tn{p), for p = 0.375 — 0.85. The scatter amongst 
the values suggested by the specific heat, order parameter, and correlation function 
data yields an uncertainty estimate for T]\f{p)- Similarly, we use the specific heat, 
cell occupancy histogram, and density-density correlation function data to estimate 
the phase-separation temperature Ti(p). 

In an effort to determine the coexisting densities with greater precision, we also 
performed grand canonical ensemble (GCE) simulations. Figure 15 shows the de- 
pendence of the density on the pair activity A. For T = 0.15 the dependence is 
smooth, but for T < 0.14 there is a discontinuity at Ac(T). Hysteresis effects are 
minimal for T > 0.12, and the coexisting densities are given by the limiting values 
as we approach Ac from above and below. The GCE studies furnish an independent 
prediction for the coexistence curve in the range T = 0.12 — 0.14. 



In Figure 16 we plot our best estimates for T/v (from canonical simulations), and 
for the coexistence curve, from both the canonical and GCE studies. Evidently, the 
coexistence curve is shifted upward in the GCE simulations, whose results place the 
tricritical point (i.e., the intersection of the lambda line and the coexistence curve), 
at around T = 0.14, as opposed to about T = 0.115 in the canonical studies. Since 
the two ensembles are not required to give identical results for finite systems, there 
is no inconsistency implied. The grand canonical results, however, seem preferable 
on several grounds. First, the shape of the coexistence curve looks more reasonable 
(more like the familiar Ising/lattice gas coexistence curve). More significantly, 
perhaps, the GCE curve actually meets the lambda line, which appears to veer 
horizontally near p = 0.4, and shows no sign of ever encountering the canonical 
coexistence curve. Given the large surface effects implicit in phase coexistence in 
a small system, the GCE, which studies homogeneous phases, seems more likely 
to give reliable results in a finite size simulation. One may speculate that the 
large surface energy depresses the transition temperature in the canonical ensemble 
simulations. A more definitive answer to this question must await the results of 
larger scale simulations. 

DISCUSSION AND SUMMARY 

We have determined the phase diagram of the lattice restricted primitive model 
in extensive Monte Carlo simulations. Our simulations reproduce the principal 
features expected on the basis of general arguments and mean-field theory, that 
is, a A-line of Neel points meeting a coexistence curve at a tricritcal point. Our 




FIGURE 15. Density versus pair activity in in grand canonical ensemble simulations {L = 16) 
for temperatures (left to right) T = 0.11, 0.12, 0.13, 0.14 and 0.15. 
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FIGURE 16. Best estimates for the location of the lambda line and coexistence curve. Filled 
symbols: canonical ensemble simulations; open symbols: GCE simulations. 



best estimate for the location of the tricritical point is Tt ~ 0.14, pt ~ 0.4, while 
mean-field theory predicts Tt = 0.299 and pt = 1/3. This comparison is consistent 
with the tendency of mean-field theories to overestimate transition temperatures, 
and underestimate transition densities (i.e., to overestimate the stability of ordered 
phases) . 

During the completion of this work we learned of Panagiotopoulos and Kumar's 
simulation results for the LRPM [19], which include a study of the effects of vary- 
ing the size of the hard core. Their results for the case of site-exclusion studied 
here support our conclusions for the phase diagram, with the exception that the 
tricritical density is somewhat larger {pt ~ 0.48 in their work) and the Neel line 
appears to be nearly vertical in the vicinity of the tricritical point. 

Because our results come almost entirely from studies of a single system size 
(L = 16), it is difficult to assign meaningful uncertainties to Tt, pt, or, indeed, the 
location of the A-line in the p-T plane. The limited comparison we made with a 
system of size L = 20 (Figure 7) suggests that the true transition temperature (for 
the A-transition) may be 5 — 10% higher than that observed for L = 16. We have 
also noted a sizeable discrepancy (~ 25%) between canonical and grand canon- 
ical simulations regarding the position of the coexistence curve along the T-axis. 
Studying larger systems is therefore the ffist order of business for future simulations 
of the LRPM. Since simulations are quite slow, with many trial moves rejected at 
the low temperatures of interest, even for the relatively small systems studied here, 
we anticipate the need for non- Metropolis simulation methods utilizing histograms, 
multi-canonical sampling, and/or cluster dynamics to augment efficiency. 

If such methods can be developed, a number of issues can be explored. It is 
of interest to understand the nature of the transitions occuring along the A-line 



and at the tricritical point, to assign these their proper universahty classes. We 
may also anticipate that the interface between the dense and dilute phases under- 
goes a roughening transition at a certain temperature. Another topic worthy of 
investigation is the global phase diagram as one goes from the Coulombic system 
studied here to a lattice gas with nearest-neighbor interactions, via a family of 
models described by a potential (e.g., of Yukawa form) having a range parameter. 
The LRPM with a lattice Green's function interaction, and on lattices that do not 
support antiferromagnetic order, are, as noted above, further topics for future work. 
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